Cyclical interactions with alliance-specific heterogeneous invasion rates 
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We study a six-species Lotka-Volterra type system on different two-dimensional lattices when each species 
has two superior and two inferior partners. The invasion rates from predator sites to a randomly chosen neigh- 
boring prey's site depend on the predator-prey pair, whereby cyclic symmetries within the two three-species 
defensive alliances are conserved. Monte Carlo simulations reveal an unexpected non-monotonous dependence 
of alliance survival on the difference of alliance-specific invasion rates. This behavior is qualitatively repro- 
duced by a four-point mean-field approximation. The study addresses fundamental problems of stability for the 
competition of two defensive alliances and thus has important implications in natural and social sciences. 
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Cyclical interactions are simple yet fascinating and power- 
ful examples of evolutionary processes yj], able to provide in- 
sights into the intriguing mechanisms of Darwinian selection 
12] as well as structural complexity [3] and pre-biotic evolu- 
tion 0]. The simplest non-trivial food web describing such 
cyclical interactions is formed by three species that have re- 
lationships analogous to the well-known rock-scissors-paper 
(RSP) game, where strategies form a closed loop of domi- 
nance. Real-life examples of such interactions include the 
mating strategy of side-blotched lizards [5], overgrowths by 
marine sessile organisms |6], and competition among differ- 
ent strains of bacteriocin-producing bacteria [7]. Cycles are 
also common in the context of evolutionary game theory JH], 
where strategic complexity [T(J often leads to RSP type of 
dominance between different strategies fiUl . 

Several theoretical aspects of multi-species cyclical dom- 
inance have already been studied in detail. For example, it 
has been established that three species in a cyclic dominance 
exhibit self-organizing behavior on the spatial grid 11121 Il3ll . 
whereby similar observations can be made also for system 
that incorporate more than three species, provided their total 
number does not exceed fourteen Ill4ll . Phase transitions and 
selection have also been studied in the predator-prey models 
allowing motion throughout inhabitable vacant sites lfl5i [Trill . 
Especially reticulate six-species models with mutation 111711 
and local mixing 111 811 have recently been studied rigorously, 
reporting the spontaneous emergence of defensive alliances 
and numerous stable spatial distributions as well as pertaining 
phase transitions in dependence on the topology of underlying 
food webs. In these systems the number of possible stationary 
states increases rapidly with the number of species because the 
solutions of subsystems (some species are missing) are also 
solutions of the whole system. The final stationary state can 
be determined by the competition of subsystem solutions (de- 
fensive alliances) that are characterized by their composition 
and spatio-temporal distribution of species. Consequently, the 
understanding of systems with many species requires the sys- 
tematic analysis of all their subsystems. In the present work 
we will study a predator-prey model that can be considered as 
a six-species subsystem of strains of bacteria using two types 
of toxins and anti-toxins in their warfare ]7[] . 

Besides increasing the number of species, the complexity 
of models can be enhanced also by the introduction of hetero- 




FTG. 1: Food web of the studied predator-prey model. Arrows point 
from predators towards prey with heterogeneous invasion rates spec- 
ified along the edges. 



geneous invasion rates between interacting individuals. Dif- 
ferences in invasion rates might affect the proportions of par- 
ticipating species in the habitat 111911 as well as geometrical 
features of patterns on the spatial grid j20ll . Moreover, it has 
recently been discovered that multi-species models of cycli- 
cal interactions comprising an even number of individuals are 
very sensitive to the independent variation of invasion rates 
12111 . The introduction of heterogeneous invasion rates also 
raises questions about the relevance of defensive alliances. 

Presently, we thus study a reticulate six-species predator- 
prey model where each site i of the square lattice is occupied 
by an individual belonging to one of the six species. Their 
corresponding distribution is given by a set of site variables 
Si = 0, . . . , 5. The predator-prey relations and the corre- 
sponding invasion rates (0 < a,(3,~f,6 < 1) are defined by 
the food web presented in Fig.Q] For this choice of parameter- 
ization the two subsystems consisting of odd and even labeled 
species are equivalent to the thoroughly studied rock-scissors- 
paper game and the system remains unchanged under cyclic 
permutation (s — > s + 2 modulo 6) of species. 

In case of homogeneous rates {a = (3 = ^ = 6 = 1) 
the system has two equivalent three-species states [denoted 
by (0 + 2 + 4) and (1 + 3 + 5)] exhibiting a self-organizing 
pattern maintained by cyclic invasions. These state are called 
defensive alliances because their members protect each other 
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cyclically against the external invaders U3,[Hl- Consider for 
example the case when the species invades allies (1 + 3 + 5), 
in particular by attacking species 1. Its intention is immedi- 
ately disabled by the species 5 that is superior to both and 1. 
Thus, the intruder is quickly abolished from the (1 + 3 + 5) 
domain by the very same species 5 that dominates species 1 
within the alliance. The same reasoning applies for all other 
possible attempts of non-allied species to invade a defensive 
alliance. Importantly, in this mechanism the proper spatio- 
temporal distribution of species plays a crucial role (the mean- 
field approximation cannot reproduce this feature). 

By introducing alliance-specific heterogeneities in the in- 
vasion rates, we are capable of analyzing the competition be- 
tween the defensive alliances. For example, we can study 
what happens when one of the associations is more aggres- 
sive towards the other (a + 1 (3) or when the internal mech- 
anism fails to assure flawless protection against the invaders 
(7 ^ S). In order to address these two issues systematically, 
it appears reasonable to introduce two parameters that, due 
to symmetries in the food web, uniquely determine the sta- 
tionary state of the system. Particularly, let G = j3 — a and 
H = j — 5 where H, G G [—1,1]. Note that the system behav- 
ior becomes trivial in two quadrants of the H — G parameter 
space because the favored defensive alliance is supported by 
both mentioned mechanisms, thus leading to its undisputed 
dominance on the spatial grid. If G < (/3 < a) and H > 
(7 > 6) allies (0 + 2 + 4) receive a two-fold advantage. First, 
due to P < a (0 + 2 + 4) invade non-members faster than in- 
dividuals in the competing alliance. Second, due to 7 > S the 
internal invasions within (0 + 2 + 4) are faster than those in 
(1 + 3 + 5), whereby faster internal cyclic invasions uphold a 
more effective protection shield against the external invaders 
by assuring a prompt response to a potential attack. An in- 
teresting competition emerges only if G > and H > (or 
equivalently if G < and H < 0), which we are going to 
explore next. Due to the symmetry of the problem our anal- 
ysis will be constrained to the parameter space spanning over 
H,G e[0,l]. 

We perform Monte Carlo (MC) simulations of the intro- 
duced six-species cyclical interaction model on the L x L 
spatial grid. Initially the six species are randomly distributed. 
The elementary steps are the following. First, two nearest 
neighbors are chosen at random, and second, if the two neigh- 
boring species form a predator-prey pair (species connected 
by an arrow in Fig. [TJ the prey is killed with the rate spec- 
ified along the arrow and an offspring of the predator occu- 
pies the prey's site. On the other hand, if the two randomly 
chosen species form a neutral pair (species not connected by 
an arrow in Fig. [TJ, or if both are identical, the second step 
dictates no action (nothing happens) and the MC simulation 
proceeds with executing step one. In accordance with the ran- 
dom sequential update, each individual is selected once on 
average during a particular Monte Carlo step (MCS). In or- 
der to characterize the stationary state we define the order pa- 
rameter m = pi + p 3 + p 5 — p — p 2 — P4, whereby p s 
(s = 0, . . . , 5) denotes the fraction of species ,s on the spatial 
grid. Here m = 1 corresponds to the complete dominance of 
the alliance (1 + 3 + 5). On the other hand, m = — 1 indicates 
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FIG. 2: Phase separation line dividing the two pure phases, charac- 
terized either by the exclusive dominance of the alliance (0 + 2 + 4) 
or (1 + 3 + 5), in dependence on H and G (/3 = 7 = 1). Linked 
squares were obtained by MC simulations while the dashed line re- 
sults from the four-point cluster mean-field approximation. 



the absolute authority of the association (0 + 2 + 4). 

We start by setting G and H equal to zero, whereby both 
alliances have an equal chance of eventually dominating the 
spatial grid via a domain growing process. In accordance with 
the above discussion of the studied model, it is reasonable to 
expect that as soon as G rises above zero (keeping H = 0) al- 
lies (1 + 3 + 5) are favored as their members invade non-allied 
species more successfully (a < (J). Thus, m = 1 in the sta- 
tionary state, meaning that the system evolves into the three- 
species self-organizing phase (1 + 3 + 5) as soon as members 
of (0 + 2 + 4) die out. However, the advantage of (1 + 3 + 5) 
given by G > can be compensated by choosing sufficiently 
large values of H, as shown in Fig. [2] As argued above, if H 
rises above zero the internal invasions within allies (1 + 3 + 5) 
slow down in comparison to (0 + 2 + 4), thus decreasing the 
effectiveness of the protection of the odd alliance and in turn 
nullifying its advantage given by G > 0. 

Strikingly though, results of MC simulations shown in 
Fig.|2]reveal a non-monotonous phase diagram in dependence 
on H. In particular, the advantage of allies (1 + 3 + 5) again 
increases if H approaches 1 [the internal invasion rate (5) van- 
ishes]. As S — > allies (1 + 3 + 5) essentially stop to invade 
each other within the alliance. In this case (S = 0) an external 
species (e.g. 2) can survive in the bulk of its neutral pair (e.g. 
in the domain of species 5). As a result, the order parameter m 
remains below 1 forming a frozen state after achieving domi- 
nance. The mechanism behind this interesting phenomenon is 
subtle and cannot be grasped at a glance. The above analysis 
was made also for some other fixed values of (3 and 7 yielding 
essentially identical results. 

Evidently, the mean-field analysis can be applied as an ana- 
lytical tool to study the behavior of the proposed predator-prey 
model. Unfortunately, the resulting master equations fail to 
confirm above results. Namely, the solution for the fractions 
of species predicts a total dominance of allies (1 + 3 + 5) 
if G > irrespective of H. This shortage of the classical 
mean-field approximation may be eliminated by applying the 
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extended versions of dynamical mean-field theory that proved 
to be very appropriate for obtaining qualitatively correct phase 
diagrams for several non-equilibrium systems. The improved 
approach involves finding a hierarchy of evolution equations 
for the configurational probabilities on fc-site clusters, where k 
characterizes the level of approximation (for details, see e.g. 
ll22l I23I0 . Nonetheless, the application of the method at the 
two-point level is still unable to account for the incursion of 
allies (0 + 2 + 4) into the G > region if H > 0. The 
first level that supports our above conjectures is the four-point 
level, of which the solution is displayed in Fig. [2] The fact that 
results of MC simulations differ from the four-point cluster 
mean-field approximation suggests that the unexpected non- 
monotonous dependence is heavily routed in the short-range 
correlation of spatial distribution, which cannot be captured 
adequately by the ansatz of a four-point level approximation. 

To obtain a better understanding of the spatial dynamics 
behind the phase diagram in Fig.|2]we examine the nature of 
phase transitions between m = — 1 and m = 1. Although 
the two competing effects of parameters G and H suggest 
that indeed a fine-tuning towards a stable co-existence of both 
alliances may be possible, thus assuring a rich diversity of 
species on the spatial grid, the reality is quite different. In 
fact, results in Fig. [2] suggest that phase transitions between 
the two pure phases are extremely sharp, and thus a stable 
co-existence of the two defensive alliances on the spatial grid 
is not feasible. Although we have only numerical arguments 
the extensive calculations performed in order to come to this 
conclusion leave extremely little room for alternatives. 

Finally, it is instructive to examine the temporal evolution 
of the order parameter in the close vicinity of the phase bound- 
ary. Figure |3 illustrates that, at the beginning, the fraction of 
even labeled species decreases on both sides of the phase sep- 
aration line. However, while above the phase boundary mem- 
bers of the even alliance go extinct, below the phase boundary 
allies (0 + 2 + 4) are able to fully recover although only a 
minute portion (as little as 0.005) of the spatial grid is occu- 
pied by its members, eventually yielding m = —1. We argue 
that this remarkable and unusual behavior, which is intimately 
linked also with the non-monotonous phase diagram presented 
in Fig. |2] is related to the formation of a suitable boundary 
layer that modifies the interaction between both defensive al- 
liances and ultimately tosses the dominance in favor of either 
(1 + 3 + 5) or (0 + 2 + 4). It is remarkable that a similar 
temporal evolution was reported in a group selection model 
rt24ll where the altruistic species almost reach extinction be- 
fore taking over the population. However, as described above, 
in our case the spatiality is a fundamental ingredient since the 
recovery of the even alliance, occurring due to its dynamical 
benefit, is an interface driven domain growing process. 

To confirm the importance of the boundary layer between 
both defensive alliances, we perform a stability analysis of 
the interface via MC simulations of the system with specially 
prepared initial conditions. In particular, we start the simula- 
tions with sharp boundaries that separate regions of initially 
randomly distributed even and odd labeled species. Within 
a 3200 x 3200 spatial grid we were able to set up 64 such 
straight boundary layers (without them interfering with each 
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FIG. 3: Temporal evolution of the order parameter for H — 0.2 
and different values of G. Solid line depicts the temporal evolu- 
tion obtained for G = 0.116, while the dashed line corresponds to 
G — 0.114. Inset shows the evolution of the order parameter from 
prepared initial states (see text for details). Lines in the inset corre- 
spond to G = 0.130, 0.117, 0.113 and 0.110 from top to bottom. 



other) and monitored in which direction the interface moved 
in dependence on G placing the system below or above the 
phase boundary shown in Fig. [2] Additionally, results were 
averaged over 10 consecutive runs to minimize unwanted fluc- 
tuations. The inset in Fig.[3]shows the results for two values of 
G below and two above the phase separation line. Note that 
to measures the difference of areas on both sides of the ini- 
tially sharp boundary between even and odd labeled species, 
thus uniquely determining also its spreading direction. Evi- 
dently, if G is set above the phase separation line in Fig. [2] the 
boundary layer moves towards the area of (0 + 2 + 4), thus 
foretelling an imminent dominance of ( 1 + 3 + 5 ) . Conversely, 
if G is set below the phase separation line the boundary layer 
moves towards (1 + 3 + 5), marking the advent of dominance 
of (0 + 2 + 4) although the total fraction of its members on the 
spatial grid at that time might be minute. Hence, if the evo- 
lution of the system is initialized from a completely random 
initial distribution of species on the spatial grid, members of 
(1 + 3 + 5) can utilize their advantage given by G > and 
rapidly start their conquest. However, as time goes by the 
seemingly defeated (0 + 2 + 4) may self-organize into small 
clusters on the spatial grid, and if G and H are set appro- 
priately, start to win back lost ground via growth along the 
interfaces, as explained above and presented in Fig. [3] 

In order to check the robustness of the above-described be- 
havior the present predator-prey model was also studied on 
the honeycomb and triangular lattices having different coordi- 
nation numbers. Figure |4] clearly shows that the honeycomb 
lattice additionally pronounces the non-monotonous variation 
of the phase boundary. Conversely, the incursion of allies 
(0 + 2 + 4) is less emphasized by the triangular lattice, which 
may be related to the increased coordination number, bringing 
the MC simulations closer to the mean-field behavior. Results 
presented in Fig. [4] further stress the importance of spatiality 
and with it related distribution of species and resulting bound- 
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FIG. 4: Phase boundaries between the stationary states (0 + 2 + 4) 
and (1 + 3 + 5) for different lattice structures. Linked symbols depict 
phase separation lines for the honeycomb (circles), square (squares), 
and triangle (triangles) lattice. 

ary layers among defensive alliances on the spatial grid. 

In sum, we have studied a six-species predator prey model 
with special heterogeneous invasion rates that were introduced 
in accordance with the two spontaneously emerging defensive 
alliances. We have shown that an increased aggression to- 



wards non-allied species could be tamed by decreasing the 
willingness of members to attack individuals within the al- 
liance itself. Remarkably though, if individuals completely 
sized to invade its own allied members the benefit of increased 
aggression towards non-allied species again increased. These 
two facts resulted in a non-monotonous dependence of al- 
liance survival on the difference of alliance-specific invasion 
rates, which we attributed to the underlying spatial dynam- 
ics of the system. We also discovered that despite the ability 
of fine-tuning two system parameters, a stable state enabling 
the co-existence of both defensive alliances on the spatial grid 
is not possible, thus resulting in sharp phase transitions be- 
tween the two absorbing states. However, the newly intro- 
duced alliance-specific heterogeneous invasion rates might yet 
prove valuable by discovering new ways of assuring biodiver- 
sity, either within generalized Lotka-Volterra models B25I1 . via 
the impact of stochasticity j26ll . or oscillatory mechanisms 
I27L We hope that the study will prove vital for the under- 
standing of the robustness of alliance formation, their com- 
petition, and for the effect of spatial structures on the evolu- 
tion of food webs in real systems, which appear to range from 
strains of bacteria to man made economic systems. 
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